In [1]:
from IPython.core.display import HTML
css_file = '../ipython_notebook_styles/ngcmstyle.css'
HTML(open(css_file, "r").read())
Out[1]:
This will be introduction:
Typical applications for GL:
Powerful method to solve PDEs over complex geometrical domains.
Very wide:
Commerical
Open source
Simulated knee flexion
Fibrous models of ligament - elastic and solid mechanics, stress models, etc.
Wrinkling analyses
Compression induced - calculation of first principal strains in, eg, paper crinkling, skin wrinkles in aging, etc. Note importance of interpolation techniques for this application.
Bi-layer structure
Importance of thickness and structure for results.
Simulation of in-plane compression of the epidermis
Another skin application; the importance of topography, swelling with moisture, etc.
Oral implants
Dentistry applications, in combination with imaging techniques.
Artificial organs and implants
Importance of mechanics, material, electronics and chemistry.
Biomedical (stents)
Combination of health and engineering aspects, with combinations of FEA and CFD.
Consumer goods
How the feel of the packaging affects the way you feel about the goods: make it feel expensive to trigger an emotional response.
EM radiation
Microwave effects on tissues, and UV-induced damage.
Movies
Physics-based computer graphics.
(Step 1 later).
Step 3:
Equivalent to a single element.
Forces $f_1$ (left end) and $f_2$ (right end).
Equilibrium equation
$$ f_1 + f_2 = 0. $$Combine with Hooke's law for the spring
$$ \begin{align} f_2 & = k (d_2 - d_1) \\ f_1 & = -k (d_2 - d_1) = f_2 \end{align} $$to get the equilibrium equations in matrix form
$$ \begin{pmatrix} f_1 \\ f_2 \end{pmatrix} = \begin{pmatrix} k & -k \\ -k & k \end{pmatrix} \begin{pmatrix} d_1 \\ d_2 \end{pmatrix} $$The stiffness matrix is symmetric but singular. The problem is that this doesn't constrain spatial translations; can take an infinite number of positions in space. Need to remember this and work around.
Go to a system with two elastic springs.
Equilibrium equation
$$ F_1 + F_2 + F_3 = 0 $$Constitutive equations
Split the original structure into elemental components
$$ {\bf f}^{(1)} = k^{(1)} {\bf d}^{(1)} $$and similarly for the second element.
Using the link between the displacements, that $d_2^{(1)} = d_1^{(2)}$, can set up a connectivity matrix linking the nodes in the different elements. Do something similarly for the forces, getting
$$ \begin{pmatrix} F_1 \\ F_2 \\ F_3 \end{pmatrix} = \begin{pmatrix} f_1^{(1)} \\ f_1^{(2)} + f_2^{(1)} \\ f_2^{(2)} \end{pmatrix} $$Expanded local equations
The location equations for each spring can be rewritten in matrix form, or they can be expanded into larger matrices and vectors. Local versions have size 2, global version has size 3.
Final version from adding up all the expanded versions
$$ K {\bf d} = {\bf f}. $$Direct assembly of the global stiffness matrix $K$.
Use the connectivity to measure the contributions or connections.
Special case: grounded spring
Partition and apply boundary condition $d_1 = 0$.
Global equilibrium equations
Q1 What is the physical meaning of $K$?
It denotes the force felt at node $i$ due to unit desplacement at node $j$ (all others fixed)
Q2 Which elements contribute to $K$?
Those between $i$ and $j$, connected to $i$.
The stiffness matrix is invertible only when suitable boundary conditions are applied.
Multiply by an arbitrary test function and integrate over the domain. Perform integration by parts to get rid of second order derivatives.
Very powerful technique - not restricted to conservative systems.
Take the strong form of the initial boundary value problem.
$$ \nabla \cdot \sigma + b = \rho \dot{v} = \rho \ddot{u}. $$Add prescribed displacement $u = \bar{u}$ and traction $t = \bar{t}$ on the boundaries, together with initial condition on displacement and velocity at $t=0$.
Consider arbitrary vector valued function $\eta = \eta(x) = \etc(\chi(X, t))$. This is the test or weighting function.
Write functional obtained by multiplying the strong form by the test function and integrating over the domain:
$$ f(u, \eta) = \int_{\Omega} ( -\nabla \cdot \sigma - b + \rho \ddot{u} ) \cdot \eta \, dv = 0. $$Expand the divergence term as
$$ \nabla \cdot \sigma \cdot \eta = \nabla (\sigma \eta) - \sigma : \nabla \eta = \nabla (\sigma \eta) - \sigma : \nabla_x \eta. $$Hence the functional becomes
$$ f(u, \eta) = \int_{\Omega} \sigma : \nabla \eta - (b - \rho \ddot{u}) \cdot \eta - \int_{\partial \Omega} \sigma \eta \cdot n \, ds $$Using that $\eta$ vanishes on the boundary you get
$$ \int_{\partial \Omega} \sigma \eta \cdot n \, ds = \int_{\partial \Omega} \bar{t} \cdot \eta \, ds. $$Finally rewrite problem as
$$ f(u, \eta) = \int_{\Omega} \left[ \sigma : \nabla \eta - (b - \rho \ddot{u})\cdot \eta \right] \, dv - \int_{\partial \Omega} \bar{t} \cdot \eta \, ds = 0, $$subject to conditions on the initial data.
Choose the test function to be a virtual displacement, $\eta = \delta u$.
The we get the principle of virtual work, and the first term in the volume integral is the internal mechanical virtual work, and the second term (less the acceleration term) and the surface integral is the external mechanical virtual work. Hence the variational principle is equivalent to the virtual work balancing off against the acceleration.
Summarizing, the weak form has the integral plus the initial conditions (in integral form).
Note that 80% of the time can be spent in meshing. In cases where the interpolation functions are NURBS or B-splines, the geometry is the mesh, so using these functions can greatly speed things up.